#!/usr/bin/perl -w
use strict;
use DBI;
use DBD::Pg;
use FindBin;
use lib ("$FindBin::Bin/..", "/net/cpp-group/Leo/bin");
use chromosomes;
use linewrap;

use db_parameters;
use ortho_species;

print <<"HEADLINE";
888888888888888888888888888888888888888888888888888888888888888888888

    table_ortho_local_paralogs

	Depends on:
        orthologs.ortholog_sets
	Updates:	
		orthologs.local_paralogs
		
888888888888888888888888888888888888888888888888888888888888888888888
HEADLINE


#
#	Connect to genome database
#
my $dbh = connect_to_panda();
use get_psql_time;
my $start_time = get_psql_time($dbh);

#_________________________________________________________________________________________
#
#	Get the relative gene positions of all orthologs for a species
#
#_________________________________________________________________________________________
sub get_ortholog_relative_gene_positions($$$\%)
{
	my ($dbh, $species, $protocol_id, $gene_id_to_pos) = @_;

	print STDERR "\t\tGet gene positions for $species:\n";
	my $sql_cmd = <<"PL/SQLCMD";
	SELECT 
        disentangled_gene_id, 
        chromosome 
    FROM 
        orthologs.ortholog_sets JOIN  
        taxon.gene_loci USING (gene_id) JOIN 
        taxon.prot_positions USING (prot_id)
    WHERE 
        protocol_id = $protocol_id AND
        species = '$species'
    ORDER BY prot_position;
PL/SQLCMD
	my $curr_chromosome = '';
	my $gene_position = -1;
	for my $data(@{$dbh->selectall_arrayref($sql_cmd)})
	{
		++$gene_position;
		if ($data->[1] ne $curr_chromosome)
		{
			$gene_position += 20000;
			$curr_chromosome = $data->[1];
		}
		$gene_id_to_pos->{$data->[0]} = $gene_position;
	}
}

#_________________________________________________________________________________________
#
#	get all ortholog sets for a species
#
#_________________________________________________________________________________________
sub get_ortholog_sets($$$\%)
{
	my ($dbh, $species, $protocol_id, $ortholog_sets) = @_;

	print STDERR  "\t\tGet ortholog sets for $species:\n";
	my $sql_cmd = <<"PL/SQLCMD";
	SELECT 
        disentangled_gene_id,
        ortholog_id
    FROM 
        orthologs.ortholog_sets
    WHERE 
        protocol_id = $protocol_id  AND
        species = '$species';
PL/SQLCMD
	for my $data(@{$dbh->selectall_arrayref($sql_cmd)})
	{
		my ($gene_id,  $ortholog_id) = @$data;
		push(@{$ortholog_sets->{$ortholog_id}},  $gene_id);
	}
}


#_________________________________________________________________________________________
#
#	get all ortholog sets for a species
#
#_________________________________________________________________________________________
sub get_disentangled_lookup($$\%)
{
	my ($dbh, $protocol_id, $disentangled) = @_;

	my $sql_cmd = <<"PL/SQLCMD";
	SELECT 
        disentangled_gene_id,
        gene_id
    FROM 
        orthologs.ortholog_sets
    WHERE 
        protocol_id = $protocol_id
PL/SQLCMD
	for my $data(@{$dbh->selectall_arrayref($sql_cmd)})
	{
		my ($disentangled_gene_id, $gene_id) = @$data;
		$disentangled->{$disentangled_gene_id} = $gene_id;
	}
}




$dbh->do("DELETE FROM orthologs.local_paralogs where protocol_id = $curr_protocol_id;");
$dbh->do("VACUUM ANALYSE orthologs.local_paralogs;");
for my $protocol (get_orthologs_protocols())
{
	my ($protocol_name,  $protocol_id) = @$protocol;
	next unless $protocol_id == $curr_protocol_id;

	my %disentangled_lookup;
	get_disentangled_lookup($dbh, $protocol_id, %disentangled_lookup);

	print STDERR "\tOrthology protocol [$protocol_name]...\n";
	for my $species ($ortho_name1, $ortho_name2)
	{
		#
		# get ortholog sets and gene positions
		#
		my %gene_id_to_pos;
		get_ortholog_relative_gene_positions($dbh, $species, $protocol_id, %gene_id_to_pos);
		my %ortholog_sets;
		get_ortholog_sets($dbh, $species, $protocol_id, %ortholog_sets);
		$dbh->do("COPY orthologs.local_paralogs(gene_id, disentangled_gene_id, ortholog_id, local_cluster_id, ".
												"protocol_id) from STDIN;") or die;
		#
		# get genes for each ortholog set and put them into local clusters depending on
		#		synteny
		#
		for my $ortholog_id (keys %ortholog_sets)
		{
			for (@{$ortholog_sets{$ortholog_id}})
			{
				die $_ unless exists $gene_id_to_pos{$_};
			}
			# get gene_ids sorted by gene position
			my @gene_ids = 	map{$_->[1]}
							sort {$a->[0] <=> $b->[0] }
							map{[$gene_id_to_pos{$_}, $_]} @{$ortholog_sets{$ortholog_id}};
			my $prev_pos = -9999;
			my $local_cluster_id = -1;
			for (@gene_ids)
			{
				my $common_gene_id = $disentangled_lookup{$_};
				die unless exists $gene_id_to_pos{$_};
				my $curr_pos = $gene_id_to_pos{$_};
				#each non-adjacent gene goes into its own cluster
				++$local_cluster_id unless ($curr_pos - $prev_pos <= 50);
				$prev_pos = $curr_pos;
				$dbh->pg_putline("$common_gene_id\t$_\t$ortholog_id\t$local_cluster_id\t$protocol_id\n");
			}
		}
		$dbh->pg_endcopy;
	}
}
	

$dbh->do ("VACUUM ANALYSE orthologs.local_paralogs");




log_pipeline_stage($dbh, $curr_protocol_id, 40, 'table local_paralogs', $start_time);
$dbh->disconnect();
print "\t",  "All finished!\n\n";




